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On the Stability of the Kuramoto Model of 
Coupled Nonlinear Oscillators^ 

Ali Jadbabaie*, Nader Motee*, and Mauricio Barahona^ 

Abstract 

We provide an analysis of the classic Kuramoto model of coupled nonlinear oscillators that goes beyond the existing 
results for all-to-all networks of identical oscillators. Our work is applicable to oscillator networks of arbitrary interconnection 
topology with uncertain natural frequencies. Using tools from spectral graph theory and control theory, we prove that for 
couplings above a critical value, the synchronized state is locally asymptotically stable, resulting in convergence of all phase 
differences to a constant value, both in the case of identical natural frequencies as well as uncertain ones. We further explain 
the behavior of the system as the number of oscillators grows to infinity. 

I. Background and introduction 

Over the past decade, considerable attention has been devoted to the problem of coordinated motion of multiple 
autonomous agents. A variety of disciplines (as diverse as ecology, the social sciences, statistical physics, computer 
graphics and, indeed, systems and control theory) are developing an understanding of how a group of moving objects 
(such as flocks of birds, schools of fish, crowds of people [11], [20], or collections of autonomous robots or unmanned 
vehicles [18], [19]) can reach a consensus and move in formation without centralized coordination. Interestingly, this has 
coincided with a surge of activity in the area of network dynamics, which focusses on the relationship between graph 
structure and dynamical behavior of large networks of diverse origin. 

A classic example of distributed coordination in physics, engineering and biology is the synchronization of arrays of 
coupled nonlinear oscillators [15], [16], [25]. Building on long-standing experiments (dating back to Huyghens and van der 
Pol), the problem of collective synchronization was explored mathematically by the Russian school of Andronov. Norbert 
Wiener [24] also recognized its ubiquity in the natural world, and even speculated about its relevance to the existence of 
characteristic rhythms in the brain [17]. 

Following on key insights by Winfree [25], Kuramoto [8] proposed in the 1970s a tractable model for oscillator 
synchronization that has become archetypal in the physics and dynamical systems literatures. (See [15] for an excellent 
review of the state-of-the-art on this model.) More recently, researchers in the control community [7], [12], [22] have 
recognized that nonlinear synchronization phenomena are mathematically related to the problem of coordination and 
consensus among multi-agent systems [6], [13]. 

II. Model description 

The classic Kuramoto model describes the dynamics of a set of N phase oscillators 6*; with natural frequencies Wj. The 
time evolution of the i-th oscillator is given by: 

K N 

ei = w 4 + -5^sin(fl i -e i ), (1) 

i=i 

where K is the coupling strength, a key parameter in the problem. One of Kuramoto's results was to show numerically 
that when the Wj's are randomly chosen from a Cauchy probability distribution in the infinite N limit, there is a critical 
value of the coupling above which all phase differences remain constant, i.e., the oscillators synchronize [8], [9]. If we 
think of the oscillators as points moving on a circle, they would rotate keeping the phase differences constant. 
Kuramoto used the magnitude R of the centroid of the points as a 'natural' measure of synchronization: 

Aerify*. (2) 

i=l 

Clearly, if all the Wj's are the same then R = 1 when all agents are in sync. If the natural frequencies are not identical 
but the oscillators synchronize, R converges to a constant R^ < 1. On the other hand, when all agents are completely 
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out of phase with respect to each other the value of R remains close to most of the time. Because it characterizes the 
dynamical behavior of the system, R is referred to as the order parameter in the physics literature. 

Kuramoto's analysis used simple trigonometry to rewrite the state equation Q in terms of the order parameter. After 
switching to a rotating frame, Eq. Q becomes: 

e i = u i -^Rwa.(fi i -il;). (3) 

In other words, each phase is modulated by the magnitude R and phase ijj of the average phasor. In physics notation, this 
constitutes a mean field or "all-to-all" model. 

With some brilliant intuition, Kuramoto showed that for an infinite number of oscillators there is a critical coupling 
K c below which the oscillators are incoherent (i.e., fully unsynchronized). In addition, there is another critical coupling 
Kl > K c above which all oscillators are synchronized. In that regime, the order parameter R grows exponentially in 
time until it saturates at a value R^K) < 1. The branch of R with K > Kl is called the fully synchronized state, while 
K < K c corresponds to the totally unsynchronized state. Kuramoto also calculated analytically the value for K c and R^ 
for a few well-known distributions in the case of an infinite number of oscillators connected all-to-all. 

Despite its success, several aspects of the well-studied TV — > oo, all-to-all Kuramoto model are still a puzzle, as 
summarized beautifully in the review by Steve Strogatz [15]. For instance, what does it mean that R stays close to zero in 
the unsynchronized state K < K c l This cannot be true at all times: when K = and the tOj's are irrational with respect to 
each other, the trajectories are dense on the TV-torus resulting in an R which will almost surely visit any number between 
and 1. However, simulations indicate that it is true most of the time. On the other extreme, the case of few oscillators 
has been tackled in the dynamical systems literature with rigorous bifurcation analysis. However, even basic results are 
not available for the large but finite TV case, which is of utmost interest in systems engineering. 

Our goal here is to perform a system theoretic analysis of the finite TV case with arbitrary connectivity. To proceed, we 
rewrite the model in terms of the incidence matrix of the undirected graph that describes the interconnection topology — 
the standard all-to-all case is then the specific case of the complete graph. We then provide several necessary as well as 
sufficient lower bounds for the critical coupling Kl- These include a bound for K below which there is no fixed-point, 
and a value of K above which there is a unique fixed-point. This extends similar results in [5] for the case of 2 oscillators 
with a finite set of values for the natural frequencies. 

III. Graph theoretical formulation of Kuramoto's model 

A good source for the necessary graph theory terminology is [4]. We formalize our results through two matrices that 
encode the topology of the connections. The incidence matrix B of an oriented graph Q a with TV vertices and e edges is 
the TV x e matrix such that: £?y = 1 if the edge j is incoming to vertex i, Bij = — 1 if edge j is outcoming from vertex 
i, and otherwise. The symmetric TV x TV matrix defined as: L = BB T is called the Laplacian of Q and is independent 
of the choice of orientation a. The Laplacian has several important properties: L is always positive semidefinite with a 
zero eigenvalue; the algebraic multiplicity of its zero eigenvalue is equal to the number of connected components in the 
graph; the TV-dimensional eigenvector associated with the zero eigenvalue is the vector of ones, ljy. It is known that the 
spectrum of the Laplacian matrix {Xi(L)} captures many topological properties of the graph. Specifically, Fiedler showed 
that the first non-zero eigenvalue MiL) (sometimes denoted the algebraic connectivity) gives a measure of connectedness 
of the graph. If we associate a positive number Wi to each edge and we form the diagonal matrix W exe := diag(Wi), 
then the matrix Lw(G) = BWB T is a weighted Laplacian which fulfills the above properties. 

In this framework, the Kuramoto model Q can be generalized to any general interconnection topology as: 

9 = w-^Bsm(B T 6), (4) 

where B is the incidence matrix of the unweighted graph, and 9 and uj are TV x 1 vectors. (It is also helpful to define 
the e x 1 vector of phase differences <fi := B T 6.) A generalization of the order parameter defined in (0 for the general 
Kuramoto model is: 

2 N 2 -2e + 21 T e cos{B T 9) 

r = w — - — - © 

It is easy to show that when the graph is complete, this is the square of the magnitude of the average phasor, i.e., for 
B = B c , we have r 2 = R 2 . While the average phasor interpretation does not generalize to general connected topologies, 
the above notion of an order parameter does generalize to arbitrary connected graphs. In fact the order parameter can be 
conveniently written in terms of the Laplacian matrix L of the underlying graph, as a measure of synchrony, or alignment. 
Specifically, after some algebra, Equation (JSJl can be written as 

r 2 = 1 - i[ e ^]*L[e J ' e ] = 1 - j- ([cos 9} T L [cos 6] + [sin 6>] T L[sin 6}). 

where [e^ ] = [e^ 1 ■ ■ ■ e : > 0N } T is the vector of complex phasors, and * denotes complex conjugate transpose. 
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The above equation provides us with an interesting interpretation of the order parameter: each individual oscillator i 
can be thought of as a rotor moving on a circle with unit radius, with velocity vector Vi = e^ 6i . The total measure of 
disagreement between all velocity vectors, or the global measure of asynchrony can be written as ^\ X^eAT W Vi ~ u ill 2 ' 
which is nothing but N 2 (l — r 2 ). In other words, the order parameter is a scaled measure of agreement among velocity 
vectors of rotors, hence a measure of synchronization. This allows us to extend the stability analysis to a graph with 
arbitrary connected topology. 

Remark 1: In the limit of small angles, the general Kuramoto model © gives the continuous-time Vicsek flocking boid 
model [21] which was analyzed in [6]: 8 = ui — (K/N) B sm(B T 8) w to — (K/N) L 8. Conversely, the classic Kuramoto 
model can be thought of as a nonlinear extension of the Vicsek model for a complete graph. 

Remark 2: It is straightforward to show that the analytical simplification (Eq.|3 in the (standard) all-to-all model appear 
as a result of the special symmetry of the Laplacian of the complete graph: 

L c = B c Bj = NI- l N lJf. (6) 



IV. Synchronization of identical coupled oscillators 

We start by considering the general Kuramoto model @ in its unperturbed version, i.e., when all the natural frequencies 
u>i are identical: 

= ~Bwa.{B T 6). (7) 

(By switching to a rotating frame, it is easily shown that we can assume that the natural frequencies u>i are all zero, 
without loss of generality.) 

Theorem 1: Consider the unperturbed Kuramoto model ([7) defined over an arbitrary connected graph with incidence 
matrix B. For any value of the coupling K > 0, all trajectories will converge to the set of equilibrium solutions. In 
particular the synchronized state is locally asymptotically stable. Moreover, the rate of approach to the synchronized state 
is no worse than (2K/irN) MiL), where M(L) is the Fiedler eigenvalue or the algebraic connectivity of the graph. 

Proof: Consider the function Ui(8) = 1 — r 2 = 4 ^ ■ s "^ jv a 2 — — , where r 2 has been defined in A simple calculation 
reveals that V e U = (2/N 2 )B s\n(B T 8) which leads to 

U(8) = \7 e U 6 = ~Yn 0T ^ < °- 

Therefore, the positive function < U(8) < 1 is a non-increasing function along the trajectories of the system. By using 
LaSalle's invariance principle we conclude that U is a Lyapunov function for the system, and that all trajectories converge 
to the set where 8 is zero, i.e., the equilibrium solutions. 

Define now the ex e diagonal matrix W(<p) := diag(sinc(0j)), where sinc(0j) = sm(<f>i) / <fii is positive for </>.; <G 
(— 7r, 7r) e . Note also that the angles move on an N- torus. Consider the proper subset of the torus in which <\> E (— tt, ir) e . 
The diagonal weight matrix W(4>) > can be thought of as phase-dependent weight functions on the graph. The trajectories 
converge to fixed-points, which are the solutions of Lw 8 := (B W((f>) i? T ) 8 = 0. The fact that (6q1n) is a locally 
exponentially stable equilibrium solution follows easily: for any connected graph the null space of the weighted Laplacian 
contains only the vector 1 at. Note that in general, for an arbitrary connected topology the system has many other equilibrium 
solutions, some of which might even be locally stable. One such example is the ring topology [10]. 

Alternatively, one could use a different Lyarrunov function similar to the approach in [14] and consider the quadratic 
Lyapunov function candidate U-2~\ 2~2i=i 2~2j=i(@i ~ ®j) 2 = 8 T L c 8, where L c = NI — 11 T is the Laplacian matrix of 
a complete graph. Note that B T 1 = 0, therefore, a simple calculation reveals that 

U 2 = -^8 T Bsm{B T 8) = -^8 T BW(<f>)B T 8 < 0. 

It is interesting to note that the Lyapunov function {/ 2 is small angle approximation of Lyapunov function U\. 

Using the same argument as above, we conclude that the largest sublevel set of {/ 2 which is contained inside \8i\ < \ l = 
1, • ■ ■ , e is positively invariant. With the quadratic function 1/2(8) however, we can show that locally, the convergence is 
exponential with the rate determined by the second smallest eigenvalue of the weighted Laplacian: 

U2 < ~\2{BW{ct>)B T )\\e^\\ 2 < -^A 2 0L)||M| 2 , 

since X 2 (BW((f>)B T ) < (2/tt)X 2 (BB t ). ■ 
Corollary 1: For the complete graph, \ 2 {L C ) = N and the synchronization rate for the mean-field model is no worse 
than 2K/-K. 

Remark 3: Similar results hold even if the topology of the graph changes in time [6]. The result can be extended to 
general notions of connectivity, i.e., when the interconnection graph is not connected at all times but there is a path 
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between any two nodes over contiguous, non-overlapping, and uniformly bounded time intervals. It is also possible to 
generalize to the case of directed graphs by introducing notions of weak connectivity [12]. 

Remark 4: The synchronization argument can be readily extended to the case of more complicated coupling functions 
/(•) (other than the sin(-) function) so long as (j> T f(cj>) > 0. 

Remark 5: The function 1^ cos(B T 6) is an energy function for the XY-model in statistical physics. It was considered 
as a Lyapunov-like function for the Kuramoto model by Van Hemmen and Wreszinski [5], as well as in [7]. 

Remark 6: The global results obtained by Watanabe and Strogatz [23] require all-to-all connectivity. An extension of 
the methodology in [23] to arbitrary topologies does not appear to be trivial. 



V. The case of non-identical oscillators 

In the rest of the paper we treat the more complicated case of oscillators with non-identical natural frequencies. Although 
there is an extensive literature for the N — ► oo case with all-to-all connectivity, we will focus here on the case of finite N 
and arbitrary topology given by Eq. l|4}. We consider the frequencies to be random perturbations which, albeit drawn from 
a probability distribution, remain constant in time, i.e., the dynamics @ is deterministic yet uncertain. This problem is 
distinct to some treatments in the physics literature, which transform the problem into a Fokker-Planck equation, effectively 
connected to a stochastic differential equation. 

Synchronization is best defined in a grounded system, where the phases are defined with respect to a reference variable 
(or 'ground'). This can be achieved by any projection Vn x (n-i) sucn that 

1 1 T L 

V T V = I, VV T = I-^j^ = ^, V T 1 N = 0. (8) 

Thus, V is a matrix of TV — 1 orthonormal vectors orthogonal to the vector ljy which generate the set of grounded 
coordinates 9 := V T 9 and frequencies uj := V t uj. The grounded Kuramoto model is: 

9 = u>- jjV T B sm{B T V6) =u>- ^V T BW{9)B T V9, (9) 

where, again, W(9) := diag(sinc(0i)) and <f> = B T V9. In this grounded system, the synchronized state is a fixed point. 

Remark 7: From Eq. l|9} it is easy to see why the natural frequencies can be centered around zero without loss of 
generality. Multiply Eq. (|9j from the left by V and use l|8} and B t 1m = to recover the original Eq. @ with new 
variables = 6 — [(to)t] In and frequencies fl = uj — (uj) In, where < uj > is the average frequency. 



VI. Bound for the asymptotic value of the order parameter 

Consider a Lyapunov function candidate based on the square of the order parameter r 2 defined in (0. The derivative 
of this function along the trajectories is 

V2 = iV 2 



— (sin B T 8) T B T B(sm B T 0) - lo t B sin B T 



which is an ellipsoid in the sin(B T 8) coordinate centered at Outside of a neighborhood of the origin given by 

\\Bsin{B T 9)\\ 2 > ^|M| 2 (10) 

the derivative is positive, resulting in growth of the order parameter. The boundary of this region contains the equilibria. 
By using an ultimate boundedness argument, the trajectories are confined to the smallest sublevel-set of r containing the 
set defined by J10I . 

We now use (II 01 to obtain an estimate of the asymptotic value of the order parameter. The vector s'm(B T 9) can be 
decomposed into two orthogonal components: y\{9), in the null space of B, and 1/2(0) in the range space of B T . The 
first component is annihilated when it is multiplied by B on the left. As a result, the region over which r 2 is positive can 
be characterized as 

Ky/\ 2 (L) 

where Aa(£) is the algebraic connectivity of the unweighted graph. We now bound the value of U over the region where 
r 2 is negative. A simple bounding reveals that 



from which 



21 T cos(B T 6») < ||1|| 2 + ||cos(B T 6»)|| 2 = 2||1|| 2 - || sin(B T e»)|| 2 , 



< jv 2 -||sini? r e|| 2 < ;v 2 HMe)l| 2 < N ' 2 - #Sl 



N 2 N 2 N 2 
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We can immediately observe that the asymptotic behavior of the order parameter is inversely proportional to the algebraic 
connectivity of the graph. Of course, because of the over-bounding, the bound is conservative — its asymptotic value is 1 
as opposed to the actual less-than-one value. Nevertheless, this gives us a bound on the growth rate of r 2 , and, as a result, 

the growth rate on r is bounded by 1 

V A 2(£) 

This means that asymptotically 



r<.l- I'" 



which would result in an increase rate of 0(— =) when the graph is complete. 

Remark 8: Consider a complete graph where the natural frequencies are independent random variables chosen from a 
normal distribution Ui ~ Af(0,a). Then ||u>||2 scales as y/Ncr, which results in a bound for r < y/l — {a /K) 2 that is 
independent of N. 

Remark 9: In [5], the authors added a linear term ui T 9 to the Lyapunov function candidate to guarantee negativity of 
the derivative everywhere except at the fixed-points, reducing the perturbed model to a gradient system. The linear term, 
however, makes the Lyapunov function indefinite. 

We will see in the next section that if K is large enough to guarantee the existence of a unique fixed point (via a 
contraction argument), condition dlOi will be trivially satisfied. This means that if K is large enough the derivative of the 
order parameter will be positive, resulting in the asymptotic stability of the synchronized state. 



VII. Bounds for the critical coupling 

As the coupling K is decreased, there is a critical value Kl below which no fixed point exists, resulting in a running 
solution for the grounded system l|9}. This means that the system cannot be fully synchronized for K < Kl- 

An easy sufficient condition for the fixed point 9* to be stable is for cf>* — B T V9* to be contained in any closed subset 
of (— §) e , which implies that \9*\ < j. This is demonstrated by taking the Jacobian of V T B smB T 9, and noting that 
it is equal to V T Bdiag[cos(B T 9*)]B T V, which is positive definite over that set. 



A. Critical value of coupling for complete graphs 

Our results generalize those of Van Hemmen et al. [5] in the case of a complete graph. Specifically, it can be shown 
that the critical value of the coupling is determined by the value of K for which the fixed point disappears. This can be 
explained by looking at the fixed point equation B sm(B T 9*) = 

Let uj max = |M|oo and note that the induced infinity norm of a matrix is the maximum absolute row sum, i.e., 
||-B||oo = dmax, where d max is the maximum degree of the graph. In the case of a complete graph, d max = N — 1. Then, 

^ _ dmax 

resulting in the following lower bound for Kl, the coupling above which a fixed point exists: 

NLOm.n:r. 



K L > 

This bound can be tightened by using the generalized inverse of V T B and bounding the component of the sm(B T 9) 
in the range of B T . The generalized inverse, denoted by (V T B)#, is equal to B T VA~ 1 , where A is the N — 1 diagonal 
matrix of the eigenvalues of the unweighted Laplacian. We therefore have the following expression 

( S m(B T 9)) R{BT) = B T VA~ 1 V T —j^r. 

Noting that L* = VA -1 V T , we have 

{Bm{B T 6)) R{B T ) = B T L#Bsm(B T 9) = B T L*^-. 

The generalized inverse of the Laplacian, in the case of a complete graph can be written as Lf = -^(1 — ^4p). Noting 
that the infinity norm of the sin vector is less than or equal to 1, and that B T L#B = B N B , we have 

||-B T w||oo < H-E^BHoo 

K - N ' 

which gives us the bound 

K L >\\B T u" 



2{N - 1) 
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This is in excellent agreement with that of Van Hemmen et al. [5] which they obtained for the simplest case of two 
oscillators. 

Remark 10: If the graph is a tree, V T B has full row rank and s'm(B9) does not have a component in the null space 
of L. In that case Kl > | \B T L#u>\ ^ is a tight bound, meaning that it is necessary and sufficient for synchronization. In 
the general case, however, this bound is just necessary. 



B. Existence and uniqueness of stable fixed points 
The fixed point equation can be written as 

(BW(B T e) B T)*^ = L*(B T e)^. 

Using Brouwer's fixed point theorem (i.e., a continuous function that maps a non-empty compact, convex set X into itself 
has at least one fixed-point), we can develop conditions which guarantee the existence (but not uniqueness) of the fixed 

' 4 ' 4 . 



point. If a fixed-point exists in any compact subset of 9 £ (— f, f ), it is stable, since this will ensure that B T 9 is between 



— § and |. We therefore have to ensure that 

K > —N max ||L&(.B T 0)||oo|M|oo. 

7T |8j|<f 

Simulations indicate that in the case of a complete graph, the infinity norm of the matrix scales as O(jj). It is 
worth mentioning that the norm of ]JL is a well studied object in the theory of Markov chains. The infinity norm of 
L^r is a measure of the sensitivity of the stationary distribution of the chain associated with L with respect to additive 
perturbations [2]. 

If the uncertain natural frequencies are 2-norm bounded, a better strategy would be to impose the boundedness condition 
with respect to the Euclidean norm. A sufficient condition for local stability of the fixed-point is for 6i to belong to 
(— j, j). This amounts to having the Euclidean norm of 9 be less than jvN. Again, using Brouwer's sufficient condition 
for existence of fixed-points we have: 

\\(BW(B T 9)B T )#\\2^^ <?ViV. 

K 4 

Hence, a sufficient condition for synchronization of all oscillators can be determined in terms of a lower bound for K: 

Rl > 4 VN\\w\\ 2 

where we used the fact that \ \{BW{B T 9)B)#\\2 = a 2 (l w ) 1 an d ^2 is the algebraic connectivity of the (weighted) graph. 
A lower bound on the minimum value of A2 occurs for the minimum value of the weight which is — . As a result, 



^> 2 «. (id 

Remark 11: Using the upper bound provided for the order parameter earlier, we can derive an upper bound for the 
asymptotic value of r at Kl- t^Kl) < Furthermore, if the stable fixed-point is in (— 7r/4, -k/£) n , then the order 
parameter is lower bounded by \/16 — 7r 2 / 4. 

C. Bounds for the existence of a unique fixed-point 

In order to guarantee the existence of a unique fixed point we use Banach's contraction principle and ensure that the 
right hand side is a contraction. By noting that the Lipschitz constant for the sinc(-) function is a s = \, we provide a 
sufficient condition for contractivity (and therefore uniqueness of the fixed-point). 

We impose the contractivity condition on the N — 1 dimensional grounded system. In the grounded case, we have 
9 = V T 9, and 

9 = (V T BW(B T 9)B T V)- 1 ^^. 
After some algebra, the contraction requirement amounts to 

^ 7T 2 NX max (L)\\w\\ 2 

K -~ W ' (12) 

where X, nax is the largest eigenvalue of the Laplacian of the graph. 

Interestingly, this value of K also ensures that the derivative of r 2 is increasing, i.e., inequality dlOl is satisfied, which 
means that the order parameter is increasing. Of course this is probably stronger than what is necessary for uniqueness, 
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as the contraction argument is only sufficient. Nevertheless, we see that there is a large enough but finite value of the 
coupling which guarantees the existence and uniqueness of fixed points. 
We now state the following theorem: 

Theorem 2: Consider the Kuramoto model for non-identical coupled oscillators with different natural frequencies Wj. 
For K > K L := 2 ^™^ , there exist at least one fixed-point for \6i\ < f or \{B T 9) l \ < f . Moreover, for K > 
tt_ n \ m ^(L)\\w\\2 jjjgjg j s on jy one g^^jg fixed-point (modulo a vector in the span of In), and the order parameter is 
strictly increasing. 

Proof: Proof of Theorem 3 : As we see before fixed point equation can be written as 

V T BW(6)B T V6 = -co (13) 



or 



9 = {V T BW{9)B T Vy l —Q (14) 
= L lv (9)- 1 ^Cu (15) 

We will use Banach's contraction principle to show that Jl 5i has a unique solution when VQ is in any compact subset 
( — § ' § ) W - ^ e tnere f° re need to show that 

N 



^-(LwihT 1 - Lwihy 1 )^ < 0*110! - 2 || (16) 
holds for some < a < 1 and some norm. Using the 2-norm, we have 

A A 

- le x it (0))x (t ce)\ wTB{w{ ° l] ~ w ^) BTv MWh 

^ JC\ (T 77 ( nJ \ VTB M W &) ~ W(9 2 )\\oo\\B T V\\ 2 \\cj\\ 2 

N 1 1 - - 

^ JF\~ Ft 7? ^ W£)IWi) - W(9 2 )\U\lj\\ 2 

^ 7F1 n <a^ \ 77 m \\ W^I|£ T IU|fli - felUlMla (17) 

a s is the Lipschitz constant of sinc(.) which is 0.5. For a weighted graph with non-negative weights W = [wij], the 
second smallest non-zero eigenvalue is given by [3] 



£i~j(/«-/j) w 



X min (L w (S))=\ 2 (L w (9)) = (N-l) inf — '■>' 13 ( i 8 ) 
with di = £^ . and D denotes the diagonal matrix with the (i, i)-th entry having value di. Since B T 6 £ (— f , f ) e 



,we have - < w^, < 1. Therefore, 



Applying l fT9l to dTTt , we get 



\2(L W (0)) > -X 2 (L) (19) 

7T 



H^i^^i)- 1 - Lw(0 2 )->||2 < f^^ll^i - fclWMIa (20) 

From here, we can find value of A"s that make this mapping contractive as follows 

N n 2 X N (L) 

k { 2 ] \ 2 JLY Mh<l (21) 



or 

^>(^)" 1 V77T 2 Z ll w ll2 (22) 



,tt 2 NX N (L) ] 

r a 2 (l) 2 



As a result for K > 2f c , where K c = (-| ) 2 ^ A ^ff IMI2, the fixed-point equation dl5> has a unique and stable solution. 



-21 A 2 (L) 
□ 
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VIII. Concluding remarks 

In this paper we provided a stability analysis for the Kuramoto model of coupled nonlinear oscillators for arbitrary 
topology. We showed that when the oscillators are identical, there are at least two Lyapunov functions which prove 
asymptotic stability of the synchronized state, when all the phase differences are bounded by We also showed that 
when the natural frequencies are not the same, there is a critical value of the coupling below which a totally synchronized 
state does not exist. Several bounds for this critical value based on norm bounded uncertain natural frequencies were 
shown to be in excellent agreement with existing bounds in the physics literature for the case of the all-to-all graph. 

We also point out that contrary to the infinite N case, there is no partially synchronized state, i.e., for values of the 
coupling below the critical value, the system of differential equations has a running solution. Furthermore, we showed 
that there is always a large enough but finite value of the coupling which results in synchronization of oscillators and 
convergence of the angles to a unique fixed-point. Another result of this paper is that the value of the order parameter is 
not zero for the critical coupling K^. In fact, at least when the fixed-point is in the (— tv/2, tt/2) region, a rough estimate 
indicates that the value of r is bounded between ^ 1& ^^ 2 « 0.62 and Future research in this direction is needed 
to determine the bound for K when the natural frequencies are not just norm bounded quantities but uncertain numbers 
chosen from a probability distribution. Finally we mention that our value for the upper bound of the order parameter is 
actually quite close to simulations. 

Our work hints at the advantageous marriage of systems and control theory and graph theory, when studying dynamical 
systems over or on networks [1]. 
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